% fig3d.m
clear; clc; close all;

%% 1. 参数设置
mu = 0.5;
x0 = 1;
q0 = 1;
k_start = 1.9;
k_end   = 2.4;
dk      = 0.002;
k_vals  = k_start : dk : k_end;

transient_steps = 300;
record_steps    = 100;

LE1_array = zeros(size(k_vals));
LE2_array = zeros(size(k_vals));

x_bif = [];
q_bif = [];
k_bif = [];
index = 1;
for k = k_vals
    x = x0;
    q = q0;
    for i = 1:transient_steps
        [x, q] = SMM(x, q, mu, k);
    end
    
    x_record = zeros(1, record_steps);
    q_record = zeros(1, record_steps);
    x_record(1) = x;
    q_record(1) = q;
    for i = 2:record_steps
        [x_record(i), q_record(i)] = SMM(x_record(i-1), q_record(i-1), mu, k);
    end
    
    x_bif = [x_bif, x_record];
    q_bif = [q_bif, q_record];
    k_bif = [k_bif, k * ones(1, record_steps)];
    
    [LE1, LE2] = LEs_SMM(mu, k, x0, q0);
    LE1_array(index) = LE1;
    LE2_array(index) = LE2;
    index = index + 1;
end

%% 绘图
figure('Name','Fig3(d) - 2D-SMM');

% 上子图：LEs 曲线
subplot(2,1,1);
hold on; grid on;
plot(k_vals, LE1_array, 'r', 'LineWidth',1.5, 'DisplayName','LE_1');
plot(k_vals, LE2_array, 'b', 'LineWidth',1.5, 'DisplayName','LE_2');
yline(0, 'k--','LineWidth',1);
xlabel('k'); ylabel('LEs');
title(['2D-SMM: \mu=', num2str(mu), ', (x_0,q_0)=(', num2str(x0), ',', num2str(q0),')']);
legend('Location','Best');

% 下子图：分岔图
subplot(2,1,2);
hold on; grid on;
scatter(k_bif, x_bif, 1, 'g', 'filled', 'DisplayName','x');
scatter(k_bif, q_bif, 1, 'm', 'filled', 'DisplayName','q');
xlabel('k'); ylabel('x, q');
legend('Location','Best');

%% 子函数
function [x_next, q_next] = SMM(x, q, mu, k)
    % 2D-SMM 映射:
    % x_{n+1} = mu*x + k*cos(q)*x
    % q_{n+1} = q + x
    x_next = mu*x + k*cos(q)*x;
    q_next = q + x;
end

function [LE1, LE2] = LEs_SMM(mu, k, x0, q0)
    MaxIter = 1000;
    Transient = 200;
    x = x0;
    q = q0;
    e1 = [1; 0];
    e2 = [0; 1];
    sum_log1 = 0;
    sum_log2 = 0;
    
    for i = 1:MaxIter
        dfdx = mu + k*cos(q);
        dfdq = -k*sin(q)*x;
        J = [ dfdx, dfdq;
              1,    1 ];
        e1_new = J * e1;
        e2_new = J * e2;
        [Q, R] = qr([e1_new e2_new], 0);
        e1 = Q(:,1);
        e2 = Q(:,2);
        if i > Transient
            sum_log1 = sum_log1 + log(abs(R(1,1)));
            sum_log2 = sum_log2 + log(abs(R(2,2)));
        end
        
        [x, q] = SMM(x, q, mu, k);
    end
    
    eff_iter = MaxIter - Transient;
    LE1 = sum_log1 / eff_iter;
    LE2 = sum_log2 / eff_iter;
end
